% RSTS_MHRUN0.M
%	samples from the prior density

% Taeyoung Doh
% created: 1/31/2008

%--------------------------------------------------------------------
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;

% MHRUN
% 0 for prior distribution
mhrun = 0;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------


nblock = 10;		% number of simulation blocks
nsimul = 10000;		% number of simulation draws in each block

%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;

% check directories
path_.mhrun = [path_.para 'mhrun' num2str(mhrun) '/'];

retdir = chkdir(path_.root);
retdir = chkdir(path_.para);
retdir = chkdir(path_.result);
retdir = chkdir(path_.mhrun);


% parameter names
pnames = feval(fnames_.pnames,msel);
npara  = size(pnames,1);


%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);
global psel mprior

%--------------------------------------------------------------------
% prior draws
%--------------------------------------------------------------------

% set clock
t0 = now;

% initialize variables
indt = 0; 
tsim = nblock*nsimul;

% reset diary file: mhrun*.log
diaryfile = [path_.result 'mhrun' num2str(mhrun) '.log'];
if exist(diaryfile,'file')
	delete(diaryfile);
end
diary(diaryfile);
diary off;

% print process
disp(' ');
disp(' ');
disp('Processing prior draws....');
disp('--------------------------');
disp(' ');
  
% loop to generate parameter draws
for iblock=1:nblock
   
    % initialize block output
	parasim = zeros(nsimul,npara);
	retcode = zeros(nsimul,1);
   
	for i=1:npara

		if prior_.mask(i)==1
			parasim(:,i) = ones(nsimul,1)*prior_.fix(i);
        else

			% beta prior
			if prior_.shape(i)==1
				a = (1-prior_.para1(i))*prior_.para1(i)^2/prior_.para2(i)^2 - prior_.para1(i);
				b = a*(1/prior_.para1(i) - 1);
				parasim(:,i) = betarnd(a,b,nsimul,1);
		
			% gamma prior		
			elseif prior_.shape(i) == 2
				b = prior_.para2(i)^2/prior_.para1(i);
				a = prior_.para1(i)/b;
				parasim(:,i) = gamrnd(a,b,nsimul,1);
		
			% gauusian prior
			elseif prior_.shape(i)==3
				a = prior_.para1(i);
				b = prior_.para2(i);
				parasim(:,i) = a + b*randn(nsimul,1);
			
			% inverse gamma prior	
			elseif prior_.shape(i)==4
				a = prior_.para1(i);
				b = prior_.para2(i);
				parasim(:,i) = sqrt( b*a^2./sum( (randn(b,nsimul)).^2 )' );

			% uniform prior	
			elseif prior_.shape(i)==5
				a = prior_.para1(i);
				b = prior_.para2(i);	
				parasim(:,i) = a + (b-a)*rand(nsimul,1);
		
			% no prior, fixed
			elseif prior_.shape(i)==0
				a = prior_.para1(i);
				parasim(:,i) = a*ones(nsimul,1);

			end
		end
	end	% i loop

   	% check determinacy of prior draws
	for j=1:nsimul
		paraj = parasim(j,:)';
if mprior==0 
alpha = paraj(1);
gamma = paraj(2);
beta  = paraj(3);
kappa = paraj(4);
tau   = paraj(5);
rho_a  = paraj(6);
rho_f  = paraj(7);
rho_i  = paraj(8);
sd_a   = paraj(9);
sd_f   = paraj(10);
sd_i   = paraj(11);
gy     = paraj(12);
lnA_0  = paraj(13);
yst    = paraj(14);

elseif mprior==1 
alpha1 = paraj(1);
alpha2 = paraj(2);
gamma1 = paraj(3);
gamma2 = paraj(4);
beta1  = paraj(5);
kappa1 = paraj(6);
tau1   = paraj(7);
rho_a  = paraj(8);
rho_f  = paraj(9);
rho_i  = paraj(10);
sd_a   = paraj(11);
sd_f   = paraj(12);
sd_i   = paraj(13);
piess  = paraj(14);
gy     = paraj(15);
lnA_0  = paraj(16);
yst    = paraj(17);
p11    = paraj(18);
p22    = paraj(19);

elseif mprior==2 

alpha = paraj(1);
gamma = paraj(2);
beta  = paraj(3);
kappa = paraj(4);
tau   = paraj(5);
rho_a  = paraj(6);
rho_f  = paraj(7);
rho_i  = paraj(8); 
sd_a_1   = paraj(9);
sd_f_1   = paraj(10);
sd_i_1   = paraj(11);
sd_a_2 = paraj(12);
sd_f_2 = paraj(13);
sd_i_2 = paraj(14);

elseif mprior == 4
alpha1 = paraj(1);
alpha2 = paraj(2);
gamma1 = paraj(3);
gamma2 = paraj(4);
beta   = paraj(5);
kappa  = paraj(6);
tau    = paraj(7);
rho_a = paraj(8);
rho_f = paraj(9);
rho_i  = paraj(10);
sd_a1  = paraj(11);
sd_f1  = paraj(12);
sd_i1   = paraj(13);
sd_a2  = paraj(14);
sd_f2  = paraj(15);
sd_i2   = paraj(16);
piess  = paraj(17); 
gy     = paraj(18);
lnA_0  = paraj(19);
yst    = paraj(20);
p11 = paraj(21);
p22 = paraj(22);

elseif mprior == 5 
alpha1 = paraj(1);
alpha2 = paraj(2);
gamma1 = paraj(3);
gamma2 = paraj(4);
beta   = paraj(5);
kappa  = paraj(6);
tau    = paraj(7);
rho_a = paraj(8);
rho_f = paraj(9);
rho_i  = paraj(10);
sd_a1  = paraj(11);
sd_f1  = paraj(12);
sd_i1   = paraj(13);
sd_a2  = paraj(14);
sd_f2  = paraj(15);
sd_i2   = paraj(16);
gy     = paraj(17);
lnA_0  = paraj(18);
yst    = paraj(19);
p11 = paraj(20);
p22 = paraj(21);

end   

%----------------------------------------------------------------------
% check the roots of the system
%----------------------------------------------------------------------
if mprior==0 | mprior==2 
A = [1 1/tau;
     0 beta]; 

B = [1+gamma/tau alpha/gamma;
     -kappa  1]; 
     
eu = 0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  determinacy condition violated  *************')
    eu = 1;
end    


elseif mprior==1 | mprior==4 | mprior==5 
A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa 0 1 0;
     0 -kappa 0 1];
     
eu = 0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
    eu = 1;
end

end 
        
  if eu==0;
   retcode(j)=0;
  else 
   retcode(j)=1;
  end 
     % print process
		if mod(j,100)==0;
			disp(' ');
			disp(['block :  ' num2str(iblock,'%5.0f')]);
			disp(['draw  :  ' num2str(j,'%5.0f')]);
			disp(' ');
		end


	indt = indt + sum(retcode==1);

    
	%save results for each block
	sname = [ path_.mhrun 'mhdraw' num2str(iblock,'%04.0f') '.mat'];
	save(sname, 'parasim', 'retcode');
    
end   
end	% iblock loop

    % elasped time
diary on;
disp(' ');
disp(        ['Start          :   ' datestr(t0)]);
disp(        ['End            :   ' datestr(now)]);
disp(        ['Elapsed Time   :   ' datestr(now-t0,13)]);
disp(sprintf('-----------------------------------------------------'));
disp(' ');
diary off;



% elasped time
diary on;
disp(' ');
disp(        ['Start          :   ' datestr(t0)]);
disp(        ['End            :   ' datestr(now)]);
disp(        ['Elapsed Time   :   ' datestr(now-t0,13)]);
disp(sprintf('-----------------------------------------------------'));
disp(' ');
diary off;


    